
/* Model estimation: EPD-NK model with price indexation */
/* By Toshihiro OKADA */



@#define prc_index=1
///////////////////////////////////////////////////////////////
//* =0: No indexization, = 1: With indexization *// 
/////////////////////////////////// ////////////////////////////

 
@#define news=1 
@#define dshock=0
@#define cshock=1
@#define hshock=1
@#define de_shock=0
@#define ce_shock=0
@#define che_shock=1
//* de_shock =1 : demand shock (general preference), teck shock, interest shock & news shock (1-4 periods) 

@#define myshock=0
//** Not used **//

@#define shock_nm=7
//** Not used **//

@#define exp_shock=4
@#define ttb_spec=8
@#define eshock=16

@#define romer=0
//** Not used **//

////// Line specification for figures ///////
lsp =1; /* setting a line specification for figures: 1='-', 2='--', 3=':' and 4='-.' */

@#define tech =1
/////////////////////////////////////////////////////////////////////////////////
//* tech=1: endogenous A, type=0: no A                                 *// 
///////////////////////////////////////////////////////////////////////////////

@#define shkm =0
@#define shki =1
@#define shku =1
@#define shkue =1
@#define shkud =1
@#define shkuc =1
@#define shkuh =1
@#define shkum =0

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//* shkm=1: money growth shock is ON, shki=1: interest rate shock is ON,  shku=1: tech shock is ON, shkud=1: preference shock is ON, shkue =1: news shock*// 
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

@#define unrealized =0
//** Not used **//

@#if tech==1
var y_g, c_g, h_g, tiv_g, w_g, a_g, tfp_g, tfpe_g,k_g,gov_g, acont_g,gg, ggn, uf pai,y,tfp, tfpe, tfpu,tfpn, u, a, acont, c,k,iv,h,w,r,prc,ppai,m,pf,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,pfn,in,lamdan,yg,hg,fpai,gap,rd,rdn,an,nf,nfn,pinf,uc,um,uh,tiv, ynobs, mc, mcn, tq, tqn, aa, aan;

 
  
 @#if ce_shock==1
   @#if eshock==8
  varexo ei eu euc  eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue;
  @#endif
  @#if eshock==9
  varexo ei eu euc eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue;
  @#endif
  @#if eshock==10
  varexo ei eu euc eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue;
  @#endif
  @#if eshock==12
  varexo ei eu euc eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11  eue;
  @#endif
  @#if eshock==15
  varexo ei eu euc eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue;
  @#endif
  @#if eshock==16
  varexo ei eu euc eue4 eue8 eue12 eue;
  @#endif
  @#if eshock==20
  varexo ei eu euc eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue15 eue16 eue17 eue18 eue19  eue;
  @#endif
 @#endif

 @#if che_shock==1
  @#if eshock==8
  varexo ei eu euc euh  eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue;
  @#endif
  @#if eshock==9
  varexo ei eu euc euh eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue;
  @#endif
  @#if eshock==10
  varexo ei eu euc euh eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue;
  @#endif
  @#if eshock==12
  varexo ei eu euc euh eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue;
  @#endif 
   @#if eshock==14
  varexo ei eu euc euh eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue;
  @#endif
  @#if eshock==15
  varexo ei eu euc euh eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue;
  @#endif
  @#if eshock==16
  varexo ei eu euc euh ef eue4 eue8 eue12 eue etfp;
  @#endif
  @#if eshock==20
  varexo ei eu euc euh eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue15 eue16 eue17 eue18 eue19 eue;
  @#endif
 @#endif

 
@#endif

@#if tech==0
var  pai,y,c,k,iv,h,w,r,prc,ppai,m,i,lamda,u,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,in,lamdan,yg,hg,fpai,gap,uc, uh, um,yobs, ynobs;
varexo ei eu eud eue euc euh  eum me men;
@#endif


 @#if tech==1
 parameters index, tau, hb  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,  gammag, gammaa,  alfa, rhoum,rhouc,rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta;
 @#endif

@#if tech==0
parameters index, hb xxx, rhoi,xpai, xy, iss, gamma, rss, yss, kss, css,wss,mss, delta, phi, theta, rho, epsilon, x, gammag, gammaa,  alfa, rhom, rhou,rhoud,rhouc, hss, s1,  sigc, sigh, sigm, ivss, lamdass, ss, notec;
@#endif

//xe=0.05;
index=0.5;
rho_eta=1.0;
hb=0.7; % Internal habit formation param.
gamma = 0.99; % discount factor  
sigc = 1.000; % 1.0, Smet=1.613 *** Need a high value to hump shaped response of i shock
sigh =2; % 0.1 or 1(CHristiano, handbook),Smet=1.265 
sigm =  2.56; % CLE=9
delta =  0.10000000/4; %capital depreciation rate
ss=6.000000001; % the second derivatives of adjustment cost: CEE=2.48, Rebello=1.3, Christiano et al (CLMR),  Nutahara =15.1 and Fujwara=2.48 (no adjustment cost: ss=0).

xxx = 0.75; %  0.8 vartheta in the paper: CLE=0.8
xpai = 1.5; % 1.5 varpi_pai in the paper: 1.5 (CLE&Gali)
xy =0.125; % 0.125 varpi_y in the paper: 0.5/4 (gali), 0.1 (CLE) 


rhoi = 0.50; % interest rate (rule) shock: 0.5 (Gali) &0.15(Smet)
rhom = 0.0;  % money growth (rulr) shock:

rhou = 0.520; % technology shock
rhoud = 0.90; % preference shock: this needs to  high to get good results
rhouc = 0.52020; % preference shock: this needs to  high to get good results
rhouh = 0.52040; % preference shock: this needs to  high to get good results
rhoum = 0.5240; % preference shock: this needs to  high to get good results
rhouf = 0.204;

phi =4.00; % the parameter related to markup rate (4.33) *.*.* 
theta = 0.33; % capital share
rho = 0.70914850000; % price stickness: 0.75 (conventional value) & rho =0 for flexible price


epsilon = 0.1/4; % R&D success probability: Comin & Gertler (0.1 per year)
gammaa = 0.750; % 0.15 technology spillover effect (Romer type R&D spec): gamma<0 for "Shoulders of a Giant" effect & gamma>0 for "Steping Toes" effect 
alfa = 0.25; % 1/(1+alfa) = elastilicity of R&D, Bransetter (0.25). This para affect hours response to productivity shocks
psi = 1-(0.1000000/4); % survival rate of A: 0.005 per year of tech. depreciation is conventional value. psi is important para for expectation shock
tau=0.18; %steady state share of gov expenditure


model(linear);
//@#if est==1

@#if ttb_spec==1
#eta1= 1;
@#endif


@#if ttb_spec==2
#eta1= ((1+rho_eta)^(-1) );
#eta2=eta1*(rho_eta);
@#endif


@#if ttb_spec==3
#eta1= ((1+rho_eta+(rho_eta^2))^(-1) );
#eta2=eta1*(rho_eta);
#eta3=eta1*(rho_eta)^2; 
@#endif


@#if ttb_spec==4
#eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3))^(-1) );
#eta2=eta1*(rho_eta);
#eta3=eta1*(rho_eta)^2; 
#eta4=eta1*(rho_eta)^3; 
@#endif

@#if ttb_spec==6
#eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5))^(-1) );
#eta2=eta1*(rho_eta);
#eta3=eta1*(rho_eta)^2; 
#eta4=eta1*(rho_eta)^3; 
#eta5=eta1*(rho_eta)^4; 
#eta6=eta1*(rho_eta)^5;
@#endif

@#if ttb_spec==8
#eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7))^(-1) );
#eta2=eta1*(rho_eta);
#eta3=eta1*(rho_eta)^2; 
#eta4=eta1*(rho_eta)^3; 
#eta5=eta1*(rho_eta)^4; 
#eta6=eta1*(rho_eta)^5;
#eta7=eta1*(rho_eta)^6;
#eta8=eta1*(rho_eta)^7;
@#endif

#iss = (1/gamma)-1;
#rss = iss + delta;
#wss = (((phi-1)/phi)*(1-theta)*(((1-theta)/theta)^(-theta))*(rss^(-theta)))^(1/(1-theta));
#yss_hss = ((theta/(1-theta))*(wss/rss))^theta;
#kss_hss = (yss_hss)^(1/theta);
#nfss= (1-psi)/epsilon;


@#if ttb_spec==8
#g1_t8 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8)  )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
#g2_t8 = (1-g1_t8)*(yss_hss) - (delta*kss_hss)- tau*(yss_hss);
#g1 = g1_t8; 
#g2 = g2_t8;
#css_hss = g2_t8; 
#etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8) );
@#endif

@#if ttb_spec==12
#g1_t12 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8) +eta9*((1+iss)^9) +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12)  )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
#g2_t12 = (1-g1_t12)*(yss_hss) - (delta*kss_hss);
#g1 = g1_t12; 
#g2 = g2_t12;
#css_hss = g2_t12; 
#etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8)  +eta9*((1+iss)^9)  +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12) );
@#endif





#hss= (wss*(css_hss^(-sigc))*((1-hb)^(-sigc))*(1-hb*gamma))^(1/(sigh+sigc)); 
#css=(css_hss)*hss;
#yss = (((theta/(1-theta))*(wss/rss))^theta)*hss;
#kss = (yss_hss^(1/theta))*hss;

#rdss = g1*yss;
#x=((nfss)^(1+alfa))/rdss; % x is (kappa_bar*delta) in the paper

#pfss = (1/(1-psi))*etas*rdss; % pfss is the steady state value of th epresent value of monopoly profits stream 
#lamdass = (1-hb*gamma)*((1-hb)*css)^(-sigc); % lamdass = the steady state velu of CHI in the paper
#mss = (((1-hb)*css)^(sigc/sigm))*(((1+iss)/iss)^(1/sigm))*((1-hb*gamma)^(-1/sigm));
#ivss = delta*kss;
#ggss =tau*yss; % ggss=staedy state value of gg, govt expnditure

 


@#if tech==0  

#iss_noa = (1/gamma)-1;
#rss_noa = iss_noa + delta;
#wss_noa = (((phi-1)/phi)*(1-theta)*(((1-theta)/theta)^(-theta))*(rss_noa^(-theta)))^(1/(1-theta));
#yss_hss_noa = ((theta/(1-theta))*(wss_noa/rss_noa))^theta;
#kss_hss_noa = (yss_hss_noa)^(1/theta);
#css_hss_noa = yss_hss_noa - delta*kss_hss_noa;
#hss_noa = (wss_noa*(css_hss_noa^(-sigc))*((1-hb)^(-sigc))*(1-hb*gamma))^(1/(sigh+sigc));
#yss_noa = (((theta/(1-theta))*(wss_noa/rss_noa))^theta)*hss_noa;
#kss_noa = (yss_hss_noa^(1/theta))*hss_noa;
#css_noa = css_hss_noa*hss_noa;
#mss_noa = (((1-hb)*css_noa)^(sigc/sigm))*(((1+iss_noa)/iss_noa)^(1/sigm))*((1-hb*gamma)^(-1/sigm));
#ivss_noa = delta_*kss_noa;
#lamdass_noa = (1-hb_*gamma)*((1-hb)*css_noa)^(-sigc);

@#endif

 


@#if tech==1

//Non-Romer type spec //
@#if romer==0
 @#if ttb_spec==1
    rd=(1+alfa)*(eta1*nf);
    pf(+1)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  eta1*iss*(1/iss)*i + eta1*(1+iss)*(pai(+1))  );    
    a = (1-psi)*nf + psi*a(-1);
 @#endif
  @#if ttb_spec==8
 rd=(1+alfa)*(eta1*nf(-7)+eta2*nf(-6)+eta3*nf(-5)+eta4*nf(-4)+eta5*nf(-3)+eta6*nf(-2)+eta7*nf(-1)+eta8*nf);
 pf(+8)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*i(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*i(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*i(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*i(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*i(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*i(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*i(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*i + eta1*(1+iss)*(prc(+7)-prc(+8)) + eta2*((1+iss)^2)*(prc(+6)-prc(+8)) + eta3*((1+iss)^3)*(prc(+5)-prc(+8)) + eta4*((1+iss)^4)*(prc(+4)-prc(+8)) + eta5*((1+iss)^5)*(prc(+3) -prc(+8)) + eta6*((1+iss)^6)*(prc(+2) -prc(+8)) + eta7*((1+iss)^7)*(prc(+1) -prc(+8)) + eta8*((1+iss)^8)*(prc -prc(+8))  );    
 a = (1-psi)*nf(-7) + psi*a(-1);
 @#endif
@#endif



c-hb*c(-1)= (1+gamma*hb)*(c(+1)-hb*c) - gamma*hb*(c(+2)-hb*c(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss/(1+iss))*(1/iss)*i + prc - prc(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*h = w - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(c - hb*c(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - uc(+1) - ud(+1) +ud );
sigm*(m-prc) = (-1/(1+iss))*(1/iss)*i + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(c-hb*c(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) - ud) );
lamda = (1-delta)*gamma*lamda(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - (1/rss)*r(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(c(+2)-hb*c(+1)) - (1/rss)*r(+1) -uc(+2) -ud(+2) ) );
iv =  (gamma/(1+gamma))*iv(+1) + (1/(1+gamma))*iv(-1) + (1/(ss*(1+gamma)))*lamda  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(c-hb*c(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) ) ) ;


yss*y = css*c + ivss*iv + rdss*rd +ggss*gg;
//y=(1+xe)*yobs-xe*me; //*measurement errors in y & x is a constant rate of ME
(1-delta)*k(-1) +delta*iv - k = 0;

tiv= (ivss*iv + rdss*rd)/(ivss+rdss); //*total investment: R&D investment+ capital investment


//y = (1/(phi-1))*a(-1) + (1/(phi-1))*u + theta*k(-1) + (1-theta)*h;
y = (1/(phi-1))*a(-1) + u + theta*k(-1) + (1-theta)*h;
h - k(-1) - (1/rss)*r + w = 0;

@#if prc_index==0
pai = gamma*(pai(+1)) + (((1-gamma*rho*psi)*(1-rho*psi))/(rho*psi))*(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*a(-1)) + (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2));
@#endif

@#if prc_index==1
pai = ( gamma/(1+gamma*rho*psi*index) )*(pai(+1)) + ( index/(1+gamma*rho*psi*index) )*pai(-1)+ (  ((1-gamma*rho*psi)*(1-rho*psi))/( (rho*psi)*(1+gamma*rho*psi*index) )  )*(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*a(-1)) + (  gamma/( (phi-1)*(1+gamma*rho*psi*index) )   )*(a - a(-1)) - (  1/( (phi-1)*(1+gamma*rho*psi*index) )  )*(a(-1) - a(-2));
@#endif


pfss*pf - (gamma*psi)*pfss*pf(+1) + (gamma*psi)*pfss*(prc - prc(+1) + iss*(1/iss)*i) = (1/phi)*yss*y +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*w -  ((phi-1)/phi)*theta*yss*(1/rss)*r;
pai = prc - prc(-1);

gg=(1/tau)*uf +y; // * govt expenditure


///////////////////////////////////////////////////////////////
// Monetary Policy (Inerest rate rule) , gove exp shock & interest rate shock //
//////////////////////////////////////////////////////////////
iss*(1/iss)*i = ((1+iss)^(xxx))*(xxx*(iss/(1+iss))*(1/iss)*i(-1) + (1-xxx)*((xpai)*(pai)+(xy)*(yg)) + v);
v = rhoi*v(-1) + ei;
uf = rhouf*uf(-1)+ef;

// Technology & news shocks //

 @#if news==1
   @#if eshock==2
   u = rhou*u(-1) +eu + eue1(-1) + eue(-2); 
   @#endif
   @#if eshock==3
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue(-3) ; 
   @#endif
   @#if eshock==4
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue(-4); 
   @#endif
   @#if eshock==5
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue(-5); 
   @#endif
   @#if eshock==6
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue(-6); 
   @#endif
   @#if eshock==7
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue(-7); 
   @#endif
   @#if eshock==8
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue(-8); 
   @#endif
   @#if eshock==9
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue(-9); 
   @#endif
   @#if eshock==10
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue(-10); 
   @#endif
   @#if eshock==11
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue(-11); 
   @#endif 
   @#if eshock==12
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue(-12); 
   @#endif
   @#if eshock==13
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue(-13); 
   @#endif
   @#if eshock==14
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue(-14); 
   @#endif
   @#if eshock==16
   //u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   u = rhou*u(-1) +eu + eue4(-4)+eue8(-8)+eue12(-12)+eue(-16); 
  
  @#endif
   @#if eshock==20
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9)+eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue16(-16)+eue17(-17)+eue18(-18)+eue19(-19)+eue(-20); 
   @#endif
 @#endif
 
 @#if news==0
   u = rhou*u(-1) +eu ; 
 @#endif


 
////////////////////////
// Preference shocks //
//////////////////////
@#if dshock==1
ud = rhoud*ud(-1) +eud; //* Preference shock
@#endif

@#if dshock==0
ud = rhoud*ud(-1); //* Preference shock
@#endif

@#if cshock==1
uc = rhouc*uc(-1) +euc; //* Preference shock
@#endif

@#if cshock==0
uc = rhouc*uc(-1); //* Preference shock
@#endif


@#if hshock==1
uh = rhouh*uh(-1) +euh; //* Preference shock
@#endif

@#if hshock==0
uh = rhouh*uh(-1); //* Preference shock
@#endif


//um = rhoum*um(-1) +eum; //* Preference shock
um = rhoum*um(-1) ; //* Preference shock






yg=y-ynobs; //* gdp gap
hg = h-hn;

pinf= (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2));
gap=(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*a(-1));
fpai = pai(+1);
ppai = pai(-1);

mc=-u + theta*(1/rss)*r +(1-theta)*w; //*real marginal cost under sticky price
tfp=(1/(phi-1))*a(-1) + u ;

tfpe=0.99999*tfpe(-1) + tfp_g + etfp;

acont=(1/(phi-1))*a(-1);
tfpu=u;
aa=a(-1);
tq=lamda-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-c+hb*c(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-c(+1)+hb*c) + uc(+1)+ud(+1) ) );

y_g=y-y(-1); //*growth rate
c_g=c-c(-1);
h_g=h-h(-1);
tiv_g=tiv-tiv(-1);
w_g=w-w(-1);
a_g=a-a(-1);
tfp_g=tfp-tfp(-1);
tfpe_g=tfp_g + etfp;
gov_g=gg-gg(-1);
acont_g=acont-acont(-1);
k_g=k-k(-1);

//////////////////////////
//* natural variables *//
////////////////////////


//Non-Romer type spec //
@#if romer==0
 @#if ttb_spec==1
    rdn=(1+alfa)*(eta1*nfn);
    pfn(+1)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  eta1*iss*(1/iss)*in + eta1*(1+iss)*(pn-pn(+1))  );    
    an = (1-psi)*nfn + psi*an(-1);
 @#endif
 @#if ttb_spec==2
 rdn=(1+alfa)*(eta1*nfn(-1)+eta2*nfn);
 pfn(+2)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2)*(1/iss)*in(+1) +   iss*((1+iss)^1)*eta2*(1/iss)*in + eta1*(1+iss)*(pn(+1)-pn(+2)) + eta2*((1+iss)^2)*(pn-pn(+2))   );    
 an = (1-psi)*nfn(-1) + psi*an(-1);
 @#endif
 @#if ttb_spec==3
 rdn=(1+alfa)*(eta1*nfn(-2)+eta2*nfn(-1)+eta3*nfn);
 pfn(+3)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3)*(1/iss)*in(+2) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3)*(1/iss)*in(+1) + iss*((1+iss)^2)*eta3*(1/iss)*in + eta1*(1+iss)*(pn(+2)-pn(+3)) + eta2*((1+iss)^2)*(pn(+1)-pn(+3)) + eta3*((1+iss)^3)*(pn-pn(+3))  );    
 an = (1-psi)*nfn(-2) + psi*an(-1);
 @#endif
 @#if ttb_spec==4
 rdn=(1+alfa)*(eta1*nfn(-3)+eta2*nfn(-2)+eta3*nfn(-1)+eta4*nfn);
 pfn(+4)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4)*(1/iss)*in(+3) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4)*(1/iss)*in(+2) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4)*(1/iss)*in(+1) + iss*((1+iss)^3)*eta4*(1/iss)*in + eta1*(1+iss)*(pn(+3)-pn(+4)) + eta2*((1+iss)^2)*(pn(+2)-pn(+4)) + eta3*((1+iss)^3)*(pn(+1)-pn(+4)) + eta4*((1+iss)^4)*(pn-pn(+4))  );    
 an = (1-psi)*nfn(-3) + psi*an(-1);
 @#endif
 @#if ttb_spec==5
 rdn=(1+alfa)*(eta1*nfn(-4)+eta2*nfn(-3)+eta3*nfn(-2)+eta4*nfn(-1)+eta5*nfn);
 pfn(+5)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5)*(1/iss)*in(+4) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5)*(1/iss)*in(+3) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5)*(1/iss)*in(+2) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5)*(1/iss)*in(+1) + iss*((1+iss)^4)*eta5*(1/iss)*in + eta1*(1+iss)*(pn(+4)-pn(+5)) + eta2*((1+iss)^2)*(pn(+3)-pn(+5)) + eta3*((1+iss)^3)*(pn(+2)-pn(+5)) + eta4*((1+iss)^4)*(pn(+1)-pn(+5)) + eta5*((1+iss)^5)*(pn -pn(+5))  );     
 an = (1-psi)*nfn(-4) + psi*an(-1);
 @#endif
 @#if ttb_spec==6
 rdn=(1+alfa)*(eta1*nfn(-5)+eta2*nfn(-4)+eta3*nfn(-3)+eta4*nfn(-2)+eta5*nfn(-1)+eta6*nfn);
 pfn(+6)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6)*(1/iss)*in(+5) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6)*(1/iss)*in(+4) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5+((1+iss)^3)*eta6)*(1/iss)*in(+3) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5+((1+iss)^2)*eta6)*(1/iss)*in(+2) + iss*((1+iss)^4)*(eta5+((1+iss)^1)*eta6)*(1/iss)*in(+1) + iss*((1+iss)^5)*(eta6)*(1/iss)*in + eta1*(1+iss)*(pn(+5)-pn(+6)) + eta2*((1+iss)^2)*(pn(+4)-pn(+6)) + eta3*((1+iss)^3)*(pn(+3)-pn(+6)) + eta4*((1+iss)^4)*(pn(+2)-pn(+6)) + eta5*((1+iss)^5)*(pn(+1) -pn(+6)) + eta6*((1+iss)^6)*(pn -pn(+6))  );     
 an = (1-psi)*nfn(-5) + psi*an(-1);
 @#endif
 @#if ttb_spec==7
 rdn=(1+alfa)*(eta1*nfn(-6)+eta2*nfn(-5)+eta3*nfn(-4)+eta4*nfn(-3)+eta5*nfn(-2)+eta6*nfn(-1)+eta7*nfn);
 pfn(+7)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7)*(1/iss)*in(+6) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7)*(1/iss)*in(+5) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7)*(1/iss)*in(+4) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7)*(1/iss)*in(+3) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7)*(1/iss)*in(+2) + iss*((1+iss)^5)*(eta6+ ((1+iss)^1)*eta7)*(1/iss)*in(+1) + iss*((1+iss)^6)*(eta7)*(1/iss)*in + eta1*(1+iss)*(pn(+6)-pn(+7)) + eta2*((1+iss)^2)*(pn(+5)-pn(+7)) + eta3*((1+iss)^3)*(pn(+4)-pn(+7)) + eta4*((1+iss)^4)*(pn(+3)-pn(+7)) + eta5*((1+iss)^5)*(pn(+2) -pn(+7)) + eta6*((1+iss)^6)*(pn(+1) -pn(+7)) + eta7*((1+iss)^7)*(pn -pn(+7))  );    
 an = (1-psi)*nfn(-6) + psi*an(-1);
 @#endif
 @#if ttb_spec==8
 rdn=(1+alfa)*(eta1*nfn(-7)+eta2*nfn(-6)+eta3*nfn(-5)+eta4*nfn(-4)+eta5*nfn(-3)+eta6*nfn(-2)+eta7*nfn(-1)+eta8*nfn);
 pfn(+8)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*in(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*in(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*in(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*in(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*in(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*in(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*in(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*in + eta1*(1+iss)*(pn(+7)-pn(+8)) + eta2*((1+iss)^2)*(pn(+6)-pn(+8)) + eta3*((1+iss)^3)*(pn(+5)-pn(+8)) + eta4*((1+iss)^4)*(pn(+4)-pn(+8)) + eta5*((1+iss)^5)*(pn(+3) -pn(+8)) + eta6*((1+iss)^6)*(pn(+2) -pn(+8)) + eta7*((1+iss)^7)*(pn(+1) -pn(+8)) + eta8*((1+iss)^8)*(pn -pn(+8))  );    
 an = (1-psi)*nfn(-7) + psi*an(-1);
 @#endif
 @#if ttb_spec==9
 rdn=(1+alfa)*(eta1*nfn(-8)+eta2*nfn(-7)+eta3*nfn(-6)+eta4*nfn(-5)+eta5*nfn(-4)+eta6*nfn(-3)+eta7*nfn(-2)+eta8*nfn(-1)+eta9*nfn);
 pfn(+9)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9)*(1/iss)*in(+8) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9)*(1/iss)*in(+7) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9)*(1/iss)*in(+6) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9)*(1/iss)*in(+5) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9)*(1/iss)*in(+4) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9)*(1/iss)*in(+3) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9)*(1/iss)*in(+2) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9)*(1/iss)*in(+1) + iss*((1+iss)^8)*(eta9)*(1/iss)*in + eta1*(1+iss)*(pn(+8)-pn(+9)) + eta2*((1+iss)^2)*(pn(+7)-pn(+9)) + eta3*((1+iss)^3)*(pn(+6)-pn(+9)) + eta4*((1+iss)^4)*(pn(+5)-pn(+9)) + eta5*((1+iss)^5)*(pn(+4) -pn(+9)) + eta6*((1+iss)^6)*(pn(+3) -pn(+9)) + eta7*((1+iss)^7)*(pn(+2) -pn(+9)) + eta8*((1+iss)^8)*(pn(+1) -pn(+9)) + eta9*((1+iss)^9)*(pn -pn(+9))  );     
 an = (1-psi)*nfn(-8) + psi*an(-1);
 @#endif
 @#if ttb_spec==10
 rdn=(1+alfa)*(eta1*nfn(-9)+eta2*nfn(-8)+eta3*nfn(-7)+eta4*nfn(-6)+eta5*nfn(-5)+eta6*nfn(-4)+eta7*nfn(-3)+eta8*nfn(-2)+eta9*nfn(-1)+eta10*nfn);
 pfn(+10)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10)*(1/iss)*in(+9) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10)*(1/iss)*in(+8) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10)*(1/iss)*in(+7) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10)*(1/iss)*in(+6) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10)*(1/iss)*in(+5) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10)*(1/iss)*in(+4) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10)*(1/iss)*in(+3) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10)*(1/iss)*in(+2) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10)*(1/iss)*in(+1) + iss*((1+iss)^9)*(eta10)*(1/iss)*in + eta1*(1+iss)*(pn(+9)-pn(+10)) + eta2*((1+iss)^2)*(pn(+8)-pn(+10)) + eta3*((1+iss)^3)*(pn(+7)-pn(+10)) + eta4*((1+iss)^4)*(pn(+6)-pn(+10)) + eta5*((1+iss)^5)*(pn(+5) -pn(+10)) + eta6*((1+iss)^6)*(pn(+4) -pn(+10)) + eta7*((1+iss)^7)*(pn(+3) -pn(+10)) + eta8*((1+iss)^8)*(pn(+2) -pn(+10)) + eta9*((1+iss)^9)*(pn(+1) -pn(+10)) + eta10*((1+iss)^10)*(pn -pn(+10))  );    
 an = (1-psi)*nfn(-9) + psi*an(-1);
 @#endif
 @#if ttb_spec==11
 rdn=(1+alfa)*(eta1*nfn(-10)+eta2*nfn(-9)+eta3*nfn(-8)+eta4*nfn(-7)+eta5*nfn(-6)+eta6*nfn(-5)+eta7*nfn(-4)+eta8*nfn(-3)+eta9*nfn(-2)+eta10*nfn(-1)+eta11*nfn);
 pfn(+11)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11)*(1/iss)*in(+10) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11)*(1/iss)*in(+9) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11)*(1/iss)*in(+8) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11)*(1/iss)*in(+7) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11)*(1/iss)*in(+6) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11)*(1/iss)*in(+5) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11)*(1/iss)*in(+4) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11)*(1/iss)*in(+3) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11)*(1/iss)*in(+2) + iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11)*(1/iss)*in(+1)  + iss*((1+iss)^10)*(eta11)*(1/iss)*in + eta1*(1+iss)*(pn(+10)-pn(+11)) + eta2*((1+iss)^2)*(pn(+9)-pn(+11)) + eta3*((1+iss)^3)*(pn(+8)-pn(+11)) + eta4*((1+iss)^4)*(pn(+7)-pn(+11)) + eta5*((1+iss)^5)*(pn(+6) -pn(+11)) + eta6*((1+iss)^6)*(pn(+5) -pn(+11)) + eta7*((1+iss)^7)*(pn(+4) -pn(+11)) + eta8*((1+iss)^8)*(pn(+3) -pn(+11)) + eta9*((1+iss)^9)*(pn(+2) -pn(+11)) + eta10*((1+iss)^10)*(pn(+1) -pn(+11)) + eta11*((1+iss)^11)*(pn -pn(+11))  );    
 an = (1-psi)*nfn(-10) + psi*an(-1);
 @#endif
 @#if ttb_spec==12
 rdn=(1+alfa)*(eta1*nfn(-11)+eta2*nfn(-10)+eta3*nfn(-9)+eta4*nfn(-8)+eta5*nfn(-7)+eta6*nfn(-6)+eta7*nfn(-5)+eta8*nfn(-4)+eta9*nfn(-3)+eta10*nfn(-2)+eta11*nfn(-1)+eta12*nfn);
 pfn(+12)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12)*(1/iss)*in(+11) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12)*(1/iss)*in(+10) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12)*(1/iss)*in(+9) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12)*(1/iss)*in(+8) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12)*(1/iss)*in(+7) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12)*(1/iss)*in(+6) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12)*(1/iss)*in(+5) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12)*(1/iss)*in(+4) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12)*(1/iss)*in(+3) + iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12)*(1/iss)*in(+2)  + iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12)*(1/iss)*in(+1)  + iss*((1+iss)^11)*(eta12)*(1/iss)*in + eta1*(1+iss)*(pn(+11)-pn(+12)) + eta2*((1+iss)^2)*(pn(+10)-pn(+12)) + eta3*((1+iss)^3)*(pn(+9)-pn(+12)) + eta4*((1+iss)^4)*(pn(+8)-pn(+12)) + eta5*((1+iss)^5)*(pn(+7) -pn(+12)) + eta6*((1+iss)^6)*(pn(+6) -pn(+12)) + eta7*((1+iss)^7)*(pn(+5) -pn(+12)) + eta8*((1+iss)^8)*(pn(+4) -pn(+12)) + eta9*((1+iss)^9)*(pn(+3) -pn(+12)) + eta10*((1+iss)^10)*(pn(+2) -pn(+12)) + eta11*((1+iss)^11)*(pn(+1) -pn(+12)) + eta12*((1+iss)^12)*(pn -pn(+12))  );    
 an = (1-psi)*nfn(-11) + psi*an(-1);
 @#endif
 
@#endif


cn-hb*cn(-1)= (1+gamma*hb)*(cn(+1)-hb*cn) - gamma*hb*(cn(+2)-hb*cn(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss/(1+iss))*(1/iss)*in + pn - pn(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*hn = wn - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(cn - hb*cn(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - uc(+1) - ud(+1) +ud );
//sigm*(mn-pn) = (-1/(1+iss))*(1/iss)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );
sigm*(0-pn) = (-1/(1+iss))*(1/iss)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );

lamdan = (1-delta)*gamma*lamdan(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - (1/rss)*rn(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(cn(+2)-hb*cn(+1)) - (1/rss)*rn(+1) -uc(+2) -ud(+2) ) );
ivn =  (gamma/(1+gamma))*ivn(+1) + (1/(1+gamma))*ivn(-1) + (1/(ss*(1+gamma)))*lamdan  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(cn-hb*cn(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) ) ) ;

 


yss*yn = css*cn + ivss*ivn + rdss*rdn +ggss*ggn;
(1-delta)*kn(-1) +delta*ivn - kn = 0;
ynobs=yn;

hn - kn(-1) - (1/rss)*rn + wn = 0;
//yn = (1/(phi-1))*an(-1) +(1/(phi-1))*u + theta*kn(-1) + (1-theta)*hn;
yn = (1/(phi-1))*an(-1) + u + theta*kn(-1) + (1-theta)*hn;

0 = (((1-gamma*0*psi)*(1-0*psi)))*(-u + theta*(1/rss)*rn +(1-theta)*wn - (1/(phi-1))*an(-1)) ;
pfss*pfn - (gamma*psi)*pfss*pfn(+1) + (gamma*psi)*pfss*(pn - pn(+1) + iss*(1/iss)*in) = (1/phi)*yss*yn +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*wn -  ((phi-1)/phi)*theta*yss*(1/rss)*rn;

pain = pn - pn(-1);
//mn = mn(-1);
//mn = 0;

tfpn=(1/(phi-1))*an(-1) + u;
aan=an(-1);
tqn=lamdan-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-cn+hb*cn(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-cn(+1)+hb*cn) + uc(+1)+ud(+1) ) );
mcn=-u + theta*(1/rss)*rn +(1-theta)*wn;

ggn= (1/tau)*uf +yn;

end;

//varobs y_g pai i c_g w_g h_g tfp_g ;
varobs y_g pai i c_g w_g h_g tfpe_g ;


//varobs y_g pai i c_g w_g h_g ;
//varobs y_g pai i c_g w_g tfp_g ;




initval;

ynobs=0;
prc = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;

a = 0;
i=0;
lamda =0;
iv=0;
tiv=0;

v=0;


nf=0;
nfn=0;


y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;
rd = 0;
pf =0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;
rdn = 0;
pfn =0;
pn = 0; 
pain = 0;
an = 0;
in=0;
ivn=0; 
lamdan=0;

hg=0;
yg=0;

tq=0;
tqn=0;
tfp=0;
tfpu=0;
tfpn=0;

aa=0;
aan=0;
mc=0;
mcn=0;

end;

steady_state_model;

 
ynobs=0;
prc = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;

a = 0;
i=0;
lamda =0;
tiv=0;

v=0;


nf=0;
nfn=0;


y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;
rd = 0;
pf =0;
//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;
rdn = 0;
pfn =0;
pn = 0; 
pain = 0;
an = 0;
in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;


 shocks;
 //var me; stderr 0.01;
 //var mee; stderr 0.01;
 var ei; stderr 0.01;
 var eu; stderr 0.01;
 var ef; stderr 0.01;
 @#if news==1
  
 

  
   @#if eshock==8
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  @#endif 


   @#if eshock==9
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  var eue8; stderr 0.01;
   @#endif 

 @#if eshock==10
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  var eue8; stderr 0.01;
  var eue9; stderr 0.01;
  @#endif 
  

  @#if eshock==12
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  var eue8; stderr 0.01;
  var eue9; stderr 0.01;
  var eue10; stderr 0.01;
  var eue11; stderr 0.01;
  @#endif 

  @#if eshock==14
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  var eue8; stderr 0.01;
  var eue9; stderr 0.01;
  var eue10; stderr 0.01;
  var eue11; stderr 0.01;
  var eue12; stderr 0.01;
  var eue13; stderr 0.01;
  @#endif 

  @#if eshock==15
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  var eue8; stderr 0.01;
  var eue9; stderr 0.01;
  var eue10; stderr 0.01;
  var eue11; stderr 0.01;
  var eue12; stderr 0.01;
  var eue13; stderr 0.01;
  var eue14; stderr 0.01;
  @#endif 

  @#if eshock==16
    var eue4; stderr 0.01;
  
  var eue8; stderr 0.01;
  
  var eue12; stderr 0.01;
  
  var eue; stderr 0.0;
  @#endif 

  @#if eshock==20
  var eue; stderr 0.01;
  var eue1; stderr 0.01;
  var eue2; stderr 0.01;
  var eue3; stderr 0.01;
  var eue4; stderr 0.01;
  var eue5; stderr 0.01;
  var eue6; stderr 0.01;
  var eue7; stderr 0.01;
  var eue8; stderr 0.01;
  var eue9; stderr 0.01;
  var eue10; stderr 0.01;
  var eue11; stderr 0.01;
  var eue12; stderr 0.01;
  var eue13; stderr 0.01;
  var eue14; stderr 0.01;
  var eue15; stderr 0.01;
  var eue16; stderr 0.01;
  var eue17; stderr 0.01;
  var eue18; stderr 0.01;
  var eue19; stderr 0.01;
  @#endif 

 @#endif
 @#if dshock==1
 var eud; stderr 0.01;
 @#endif
 
 @#if cshock==1
 var euc; stderr 0.01;
 @#endif
 
 //var euc; stderr 0.1;
 
@#if hshock==1
 var euh; stderr 0.01;
@#endif
 //var eum; stderr 0.1;

var etfp; stderr 0.0010;

end;


estimated_params;

index,0.3291,0.01,0.99,BETA_PDF,0.5,0.25;

rho_eta,,0.0001,, normal_pdf, 1, 0.2;
xxx,,0.0001,0.9999, beta_pdf, 0.75, 0.1;
xpai,1.7985,0.0001,, normal_pdf, 1.5, 0.25;
xy,0.0893,0.0001,,normal_pdf, 0.125,0.05; 


hb,0.4205,0.0001,0.99999,BETA_PDF,0.7,0.1;
rho,0.7813,0.0001,0.99999,BETA_PDF,0.5,0.1;
//sigh,1,0.0001,,normal_pdf, 2, 0.75;
sigc,1,0.25,3,normal_pdf, 1.5, 0.375;
//sigm,3,0.0001,,normal_pdf, 2.56, 0.5;
ss,7,0,15,normal_pdf, 4, 1.5;

phi,4,1.50000,,normal_pdf, 6, 0.900;
//psi,0.975,0.0001,,beta_pdf,0.975,0.02;
psi,0.975,0.00001,0.99999,beta_pdf,0.975,0.02;

rhou,0.5,0.00001,0.99999,beta_pdf, 0.5, 0.2;
rhoi,0.5,0.00001,0.9999,beta_pdf, 0.5, 0.2;
//rhouf,0.5,0.00001,0.9999,beta_pdf, 0.5, 0.2;

@#if dshock==1
rhoud,0.5,0.00001,0.9999, beta_pdf, 0.5, 0.2;
@#endif

@#if hshock==1
rhouh,0.5,0.00001,0.99999, beta_pdf, 0.5, 0.2;
@#endif

@#if cshock==1
rhouc,0.5,0.00001,0.99999, beta_pdf, 0.5, 0.2;
@#endif 

rhouf,0.5,0.00001,0.99999,beta_pdf, 0.5, 0.2;


alfa,0.3,0.0001,, normal_pdf, 0.25, 0.1;
//alfa,0.3,0.0001,, normal_pdf, 0.25, 0.125;
//alfa,0.1,,,inv_gamma_pdf, 0.15, 4;
//alfa,0.1,,,gamma_pdf, 0.25, 2;

//stderr me,0.1,,,inv_gamma_pdf, 0.05, 0.02;
stderr eu,0.1,,,inv_gamma_pdf, 0.008, inf; 
stderr ei,0.1,,,inv_gamma_pdf, 0.003 ,inf;
//stderr ef,0.1,,,inv_gamma_pdf, 0.5, inf;
@#if cshock==1
stderr euc,0.1,,, inv_gamma_pdf, 0.01, inf;
@#endif
@#if hshock==1
stderr euh,0.1,,,inv_gamma_pdf, 0.01, inf;
@#endif
@#if dshock==1
stderr eud,0.1,,,inv_gamma_pdf, 0.01, inf;
@#endif
stderr ef,0.1,,,inv_gamma_pdf, 0.01, inf;


//stderr tfp_g,0.1,,,inv_gamma_pdf, 0.01, inf;
//stderr k_g,0.1,,,inv_gamma_pdf, 0.005, inf;


@#if news==1
 
 
@#if eshock==16
 stderr eue4,0.01,,,inv_gamma_pdf,0.004,inf;
 stderr eue8,0.01,,,inv_gamma_pdf,  0.004,inf;
 stderr eue12,0.01,,,inv_gamma_pdf,  0.004,inf;
 stderr eue,0.01,,,inv_gamma_pdf,  0.004,inf;

//stderr etfp,0.1,0,0.005,normal_pdf, 0.0008, 0.0001;
stderr etfp,0.1,0,,gamma_pdf,  0.0016, 0.0003;
//stderr etfp,0.1,0,,gamma_pdf, 0.0012, 0.0002;
//stderr etfp,0.001,,,uniform_pdf, , ,0,0.005;



corr eue4, eue8, 0.24, , , beta_pdf, 0.0, 0.45, -1, 1;
corr eue4, eue12,-0.27614, , , beta_pdf, 0.0, 0.45, -1, 1;
corr eue4, eue, -0.49, , , beta_pdf,0.0, 0.45, -1, 1;

//corr eue4, eue8, 0.5, , ,uniform_pdf, , ,-1,1;
//corr eue4, eue12, 0.37372, , ,uniform_pdf, , ,-1,1;
//corr eue4, eue, -0.17372, , ,uniform_pdf, , ,-1,1;

//corr eue4, eue12, 0.0615, , ,uniform_pdf, , ,-1,1;
//corr eue4, eue, -0.4558, , ,uniform_pdf, , ,-1,1;
//corr eue4, eu, 0.0041320, , ,uniform_pdf, , ,-1,1;

corr eue8, eue12, 0.394, , , beta_pdf,0.0, 0.45, -1, 1;
corr eue8, eue, 0.2114, , , beta_pdf,0.0, 0.45, -1, 1;

//corr eue8, eue12, 0.4378, , ,uniform_pdf, , ,-1,1;
//corr eue8, eue, 0.0704, , ,uniform_pdf, , ,-1,1;
//corr eue8, eu, -0.610320, , ,uniform_pdf, , ,-1,1;

corr eue12, eue, 0.19, , , beta_pdf, 0.0, 0.45, -1, 1;
//corr eue12, eue, 0.38138, , ,uniform_pdf, , ,-1,1;
//corr eue12, eu, 0.0041320, , ,uniform_pdf, , ,-1,1;

//corr eue, eu, -0.028320, , ,uniform_pdf, , ,-1,1;

@#endif






@#endif


end;

estimated_params_init;


index,0.082150410525124;
rho_eta,1.153143392;
xxx, 0.572872535052;//* 0.8
xpai, 0.85908788908279834698058;//*1.5596;
xy, 0.3427;//*0.1005; 
hb,0.8878792188639;//*0.5006
rho,0.73576;
sigh, 2.00;//*1.9;
sigc,0.9986990000;//*1.9904;
sigm,2.56;
ss, 6.5200;//*4.23;
phi, 3.1505479157;//*4.2783;
psi,0.9708658605755;

rhou,0.90162519357;
rhoi,0.460324640;//*0.2081;//*0.4940;
rhouc,0.0015030896058;//*0.5981;//*0.7188;
rhoud,0.5578500;//*0.5981;//*0.7188;
rhouh,0.947559;//*0.3039;
rhouf,0.94628621;//*0.3039;


alfa,0.2719520412081;


//stderr me,0.0055;
stderr eu,0.00800796312050; 
stderr ei,0.002510328;
stderr ef,0.0056510539;



stderr etfp, 0.0038234;



//stderr tfp_g,0.00831;
//stderr k_g,0.000651;

@#if dshock==1
stderr eud,0.0574200121041024;
@#endif

@#if cshock==1
stderr euc,0.0414581001215083106;
@#endif

@#if hshock==1
stderr euh,0.0165068501219050190;
@#endif

@#if news==1
 

 @#if eshock==16
  stderr eue4,  0.003001839510022;
  stderr eue8, 0.0031045180032;
 stderr eue12, 0.003092045150021;
  stderr eue,  0.002169520016;
 @#endif

 
@#endif


end;

 
 
//options_.gmhmaxlik.iterations = 10; 
//options_.gmhmaxlik.number = 10000000;
options_.gmhmaxlik.nscale = 1000000;
//options_.gmhmaxlik.nclimb = 10000;


@#if che_shock==1 
//estimation(mode_compute=6, mode_file=merr_growth1noAdj_index_new03mm2new_mode,optim=('NumberOfMh',4,'ncov-mh',70000,'nscale-mh',500000, 'AcceptanceRateTarget',0.25),filtered_vars, smoother, moments_varendo, conditional_variance_decomposition = [ 4 8 12 16 20 32 40 60], plot_priors = 0,prior_trunc=0,bayesian_irf,irf=40, irf_shocks=(eue4,eue8,eue12,eue, eu,ei,euc, euh, ef), datafile=data_growth,prefilter=0,first_obs=1,presample=4, mh_replic=300000,mh_nblocks=2,mh_drop=0.7) y_g c_g h_g w_g pai i tiv_g tfp_g tfpe_g a_g gov_g acont_g y prc a acont tfp tfpe tfpu c h tiv iv rd gg k w  r mc ;
estimation(mode_compute=0, mode_file=merr_growth1noAdj_index_new03mm2new_mode,optim=('NumberOfMh',4,'ncov-mh',70000,'nscale-mh',500000, 'AcceptanceRateTarget',0.25),filtered_vars, smoother, moments_varendo, conditional_variance_decomposition = [ 4 8 12 16 20 32 40 60], plot_priors = 0,prior_trunc=0,bayesian_irf,irf=40, irf_shocks=(eue4,eue8,eue12,eue, eu,ei,euc, euh, ef), datafile=data_growth,prefilter=0,first_obs=1,presample=4, mh_replic=300000,mh_nblocks=2,mh_drop=0.7) y_g c_g h_g w_g pai i tiv_g tfp_g tfpe_g a_g gov_g acont_g y prc a acont tfp tfpe tfpu c h tiv iv rd gg k w  r mc ;


@#endif
 

@#endif

verbatim;
options_.datafile = 'data_growth';
var_list_ = char('tfp_g', 'tfpe_g', 'a_g', 'acont_g ', 'a', 'acont', 'tfp', 'tfpe', 'tfpu');
options_.smoother = 1;
options_.order = 1;
[oo_,M_,options_,bayestopt_]=evaluate_smoother('calibration',var_list_,M_,oo_,options_,bayestopt_,estim_params_);
end;

